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We present a new approach for spatiotemporal focusing through complex scattering media by 
wave front shaping. Using a nonlinear feedback signal to shape the incident pulsed wave front, we 
show that the limit of a spatiotemporal matched filter can be achieved, i.e., the wave amplitude at 
the intended time and focus position is maximized for a given input energy. It is exactly what is also 
achieved with time-reversal. Demonstrated with ultrasound experiments, our method is generally 
■ applicable to all types of waves. 

<D ' 
(N 



o 



I. INTRODUCTION 



Most imaging systems are based on the ability to focus a wave beam inside the area of interest. In the context 
of echographic imaging, focusing of ultrasound in the human body can be achieved with a transducer array and 
electronic delay lines: the same pulsed waveform is sent from each transducer with the appropriate delay making all 
the waveforms converge in synchrony at the desired focus location. However, that principle is not practicable any 
^ | more as soon as the sample thickness becomes larger than the mean free path, i.e., the average distance between 
O ■ two random scattering events. In such a strongly scattering medium, it has been shown with ultrasonic waves that 
spatiotemporal focusing can be achieved using time-reversal: a training pulse is first sent from a source located at 
the intended focal point, travels through the scattering medium and is captured at a transducer array, the 'time 
' reversal mirror' The waveforms received on the time reversal mirror are flipped in time and sent back, resulting 
in a wave converging at the desired focus location. Time-reversal focusing is optimal is the sense that it achieves a 
spatiotemporal matched filter, a term well-known from signal processing; for a given input energy, the amplitude of 
the pulse at the focal spot is maximal [2]. Time-reversal has also been implemented for microwaves Q. For optical 
waves, time-reversal has not been demonstrated yet, since optical time-reversal mirrors are challenging to realize. 

Therefore in optics a different approach, named 'wave front shaping' (WFS), has been taken to steer light through 
strongly scattering media The principle is based on spatial phase modulation by Spatial Light Modulators (SLM), 
' which enables one to modulate the complex amplitude of the independent modes of the incident wave front. The 
intensity in the selected output mode, e.g. the focal spot, is used to match the wave front to the scattering medium 
by an adaptive algorithm. Owing to their enormous amount of pixels, state-of-the-art spatial light modulators allow 
achieving transmission of substantial amounts of the input energy [5] . 

A big advantage of adaptive wave front shaping approaches lies in the fact that a direct access to the amplitude 
• • ' at the focus is not required, such that any type of intensity probe can be used. For instance the fluorescence from 
dye molecules allows to focus on objects hidden inside a complex medium The concept has been applied to other 
types of waves, such as surface waves 0] and ultrasonic waves in the single scattering regime for medical imaging Q 
j_J \ and MR-guided ultrasonic therapy 

Initially introduced in the monochromatic regime, WFS has recently been extended to broadband light for spa- 
tiotemporal focusing [l(| • In this work it was shown that spatial control of the input waves is sufficient to achieve 
spatiotemporal control of the scattered waves. The key idea is to use coherence gating for optimizing a single specific 
predetermined point in time. Performing the optimization at a specific time in the speckle amounts to combine all 
the scattering paths that have about the same length. This principle also works efficiently without defining a focusing 
time in the case when the temporal delay spread induced by the scattering medium is comparable to the input pulse 
duration [llj. In both cases, the missing control of temporal degrees of freedom is compensated by manipulating 
spatial ones. In a complementary way, temporal control of the input light based on pulse shaping [12fl can be used 
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Figure 1: Experimental system. From N points in the emitter plane on the left, broadband signals are transmitted 
through a two-dimensional multiple scattering medium and recorded at M points in the receiver plane on the right. 

Details of the experiment are specified in section IIIII 



to achieve spatiotemporal focusing by spatially localized phase compensation [l3| . Both approaches exploit the fact 
that temporal and spatial degrees of freedom are mutually convertible in a complex medium [14j . 

In the present paper, we demonstrate how to fully control both spatial and temporal degrees of freedom to achieve 
optimal spatiotemporal focusing for a broadband ultrasonic wave propagating through a complex multiple scattering 
medium. Like in the context of time-reversal focusing, we use the term 'optimal' in the sense of a spatiotemporal 
matched filter, i.e., for a given input energy, our method maximizes the wave amplitude at the intended time and 
focus position. The difference with time reversal is that no amplitude measurement is required at the focus location. 
Instead spatial and frequency resolved wave front shaping can be achieved with a nonlinear feedback intensity signal 
and is capable to reach optimal spatiotemporal focusing. Thus our concept should be generalizable to all types of 
waves. 



II. THE SPATIOTEMPORAL MATCHED FILTER APPROACH 



A. Optimal focusing by time reversal 

The experimental system is shown in Fig. [TJ Our goal is to maximize the signal amplitude at some focus point mo 
in the receiver plane at time to for a given energy radiated from N points in the emitter plane. In signal processing, 
this is known as the spatiotemporal matched filter approach. 

a. Spatial matched filter At first we recapitulate the conditions for optimal spatial focusing in the sense of a 
spatial matched filterQ- Let (E(u))) n , 1 ^ n 5^ N be the input signals in the emitter plane, and correspondingly 
(F(uj)) m the signals received in the image plane. The propagation through the medium from the emitter points to 
the control points at one frequency is described by: 

F{u)=T(u)E(u), (1) 

where T(w) is the propagation operator, or transfer matrix, between the points at a given frequency u). The amplitude 
received at the focal point mo at a single frequency uj is given by 

F mo (E) = m = pi , (2) 

where 

= {0,..., 0,1,0,...., 0} (3) 
is the projection onto mo- The inequality of Cauchy-Schwartz sets the upper bound for this expression: 

(_E|*T*F(°)) 



< 



; T *^(0) 



(4) 



\\E\\ 

Equality holds, if E ex* T*F®\ or 

Vn E n {uj) = a(u})T* lon (u}). (5) 

The maximum amplitude at the focal point mo is reached when all E n are proportional to the complex conjugate of the 
respective transmission coefficient. This situation corresponds to a phase conjugation experiment or a time-reversal 
experiment at a single frequency. 
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b. Temporal matched filter Let us denote by T mon (t) the impulse response between emitter n and receiver mo, 
which is the Fourier Transform of the respective transfer function in the frequency domain T mon (uj). When one aims 
at maximizing the amplitude received at time to = 0, it is well-known that E n (t) — T mon (—t) has to be sent from n. 
The signal at mo then writes as 

F mo (t) = E T ^n(t) ® T mon (-t) =J2 I du e iut T mon (u;)T^ on (u;) (6) 

n n •* 

which is the autocorrelation of the impulse response, or the inverse Fourier transform of the power spectrum, re- 
spectively. It is an even function with its maximum value at t — when all frequency contributions add in phase. 
Hence the magnitude at to is directly proportional to the transmitted energy. As T rnon (—t) is the Fourier Transform 
of T* n (to), it also achieves the spatial matched filter (Eq. [5]). 

This signal exactly corresponds to the one transmitted in a broadband time-reversal experiment. In the first step, 
the control point mo in the image plane behaves like a source, such that the wave field recorded by the two-way 
transducers in the emitter plane is given by T nmo (t). In the second step, the recorded fields are emitted in a time- 
reversed manner. Under assumption of reciprocity the emission writes as E n (t) — T raTOo (— t) = T mo7l {— t) and thus the 
focus is given by Eq. [5] Hence, time-reversal exactly achieves a spatiotemporal matched filter since it implements a 
spatial filter at each frequency of the incident pulse and the correct phase relation between the frequency components 
is intrinsically recovered. 

B. Optimal focusing by wave front shaping 

In the following we lay out the wave front shaping algorithm and investigate its performance as a matched filter. In 
order to achieve optimal spatiotemporal focusing of a broadband pulse, we perform the wave front shaping algorithm 
in two stages. During the first phase, we achieve optimal spatial focusing. Based on this result, optimal temporal 
focusing is achieved after the second phase. 

1. Optimal spatial focusing 

The direct implementation of wave front shaping as applied in optics to monochromatic acoustic waves would be 
described as follows: at first an arbitrary wave front is sent from the emitter plane. The intensity is then measured at 
the intended focal point. This can be done directly in measuring the pressure field with a hydrophone or, in the case 
of a soft elastic medium, in measuringthe displacement induced at the focal point by the radiation force, which is 
proportional to the acoustic intensity [8fl . Then the phase of the first transducer is cycled between and 2n while the 
intensity is recorded. The phase value for which intensity at focus is found maximum is then stored in the memory 
and the same operation is repeated for each transducer. Once the optimal phase has been stored for each transducer, 
the wave front can be synthesized by exciting each transducer with the determined optimal phase. 

a. Hadamard method Instead of using such an optimization scheme element by element, we perform a basis 
transformation to construct virtual transmitters Ej, defined as linear combinations of real ones, which greatly improves 
the sensitivity of the method and reduces the number or iterations Q . The coefficients of the combination are taken 
as the columns H n (l < n < N) of an N by N Hadamard matrix H, with elements Hj n e {—1; 1}. This choice ensures 
that the amplitude transmitted from each transducer is maximal. The basis transformation is performed by 

E = H~ 1 E (7) 

and the inverse transformation by 

E = HE. (8) 

The field in the receiver plane is given by 

F = TE = THE = TE (9) 

where T is the transfer matrix between the emitters in the hadamard basis and the receivers in the canonical basis. 
For the optimization algorithm, we arbitraritly chose the first column as the reference. AN coded beams are sent into 
the scattering medium as: 
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Cj,s =E ref + E j e iS (10) 

je[l,N}; 5e[0,7r,|,-|] (11) 

After each transmission the resulting intensity at the focal point mo is given by 

Ij,S = a lef + a j + 2a re/ a., cos(0 re/ - <f>j - 5), (12) 

where a^e 1 ^ = T mo j is the corresponding transmission coefficient from the transfer matrix T. From the set of 
measured intensities we are now able to calculate the optimal wave front. The complex amplitude to be transmitted 
from each virtual transducer to optimize the intensity is given by the following expressions: 



i smf I j | j Iref 

? sm f — n sm S 'i4>T ) a j — a — V 1 3-° + ~~ 



R^= fl ^e'*i 7 J V^> u 1 2 (13) 

|^/ = re/ -^, = arctan(i^|^) 

where the (/>^ m ^ compensate for the phase differences between the contributions at the focal point. The choice of a™^ 
ensures that WFS achieves a spatial matched filter. Consequently, all £^ m ^are equivalent to their corresponding 
up to the global phase difference </> r e/, 

where selects the mo-th row of the transmission matrix. Finally, a change of basis leads back to the transducer 
basis in the emitter plane 

E smf = nE srnf = ^ref h(T*)*.F(°) = Ne**™* (T*)*F (0) , (15) 

showing, that the wave front shaping algorithm fulfils the condition of the spatially matched filter (Eq. [5]). The 
emission is equivalent to the time reversal experiment. With the definition of the Hadamard matrix, HH* = Nl, the 
prefactor N arises since we omitted the equivalent transformation in the receiver plane. Since this proportionality 
factor does not impair any conclusions (see Eq. [5]), we omit it in the following. 

b. Enhancement factor After optimization for focusing on point mo, the emitted signal from transducer n is 

E srnf = _ T *n n (lg) 

where the denominator normalizes the transmission. The received signal at focus point mo is 



*m = ] = - J / y \J-m n\ ■ (I'J 



ITLqU I 



The average energy received at the focus point mo is 

J (1) = (\F m f) = (j2 \ T ^on\ 2 ^ =N(\T rnon \ 2 ) , (18) 

where the brackets denote the statistical average. We assume that the transmission coefficients are independent 
random variables and follow a circular Gaussian distribution. Without optimization, the normalized initial emission 
is 

E n = , 1 = -i- (19) 

The average received signal on point mo is 

J"° = (lF„J 2 ) = ( £ " r " , °" E ": T ' ; °"' \ = (|r,„„l 2 ) . (20) 
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The average enhancement, which also gives the signal to background ratio to other non-optimized modes, is therefore 
given by 

«=j®= N - ( 21 ) 

Initial works in optics [ij used spatial light modulators which were limited to phase-only control of the emitted 
amplitudes, such that the emission is given by 

E%> = ^fpPH (22) 

V Jy \ 1 m n\ 

The amplitude at the focus is 

Fm = —7= J]] |Tmo"l • (23) 

* n 

In this case the received energy at focus point m is 

J {2) = (\F mo \ 2 ) = (\T mon \ 2 ) + (N-1) (\T mon \f (24) 
Hence the enhancement is lowered by the well-known prefactor 

(\T man \) 2 TT 



= - « 0.8. (25) 



c. Algorithm for optimal spatial focusing During the first phase, we achieve spatial focusing for all frequency 
components contained in emission spectrum of the broadband pulse with its bandwidth Aujbb- Since the transmission 
coefficients T nm (oj) (Eq. [1]) are random variables of frequency with a correlation length 5lo < AiObb, the optimal 
emission E sm $ (oj) (Eq. [TSf has to be determined frequency-resolved for intervals Auj < Suj. Consequently, the 
Hadamard algorithm is performed frequency by frequency in steps of Auj, each time using a narrowband pulse with a 
Gaussian spectrum of a bandwidth Alu for the emission. When afterwards the full broadband signal with the optimal 
coefficients E sm f (u>) is emitted, spatially matched focusing is achieved. Contrary to the time-reversal experiment, 
temporal focusing can not yet be expected. For this, all frequency components have to be in phase, which is intrinsically 
fulfilled by time reversal. However, after the wave front shaping optimization, the unknown phase factor e^™*^' 



2. Optimal temporal focusing 

During the second phase, the spatially focused signal is additionally focused in the time-domain. For that, the 
phase factor e l< >' re f^ J > needs to be determined and its conjugate equally multiplied to all emissions £* m/ (w). Such 
spatially invariant phase filtering leads only to a temporal redistribution of the signal at the focal point. The time- 
integrated intensity at the focus is unaltered by the filtering process, excluding it as an appropriate feedback signal. 
Instead of detecting the linear intensity at the focal spot, we place a nonlinear detector which allows us to measure the 
time-integrated second harmonic intensity. This detector is treated in more details in the following paragraph. The 
second harmonic energy is sensitive for a temporal redistribution of the pulse energy. Used as a feedback signal, it 
reaches its maximum when the transmitted signal is shortest in duration, corresponding to the so-called Fourier-limit, 
when the waves at all frequencies are in phases at the focus. This technique has similarly been used in optics for the 
compression of dispersed ultrashort laser pulses [l5| . 

The most straight-forward approach is an adaptive algorithm equivalent to the spatial optimization scheme described 
in the beginning of the section. The frequency range of the broadband emission spectrum is subdivided into intervals 
Auj < 5lo. Within each iteration j of the optimization, we emit the full signal as determined for optimal spatial 
focusing, E^ 71 ? , but modified with a spatial invariant phase mask e lA *^^. Step by step, all frequency intervals are 
addressed successively. Before the first iteration, the phase mask is A(f)(uj) — OVw. In step j, the phase mask is 
modified only at the frequency ujj. We perform several emissions with the phase A(f)(ujj) increasing from A<j){ujj) = 
to A<p(ujj) = 2tt, while the second harmonic energy Jj H is recorded. As a function of each A<fi(ujj), J^ H follows 
a cosine behaviour J? H oc cos[A(p(ujj) + tfisH (w)] with a phase offset 4>sH(uj), which depends on the actual phases 
of all other frequency components. The optimum phase offset A<fi(uij) — ~<fisH(uj) is found when Jj reaches its 
maximum. The modified phase mask is subsequently used in the next step address, in which the next frequency 
interval u)j+i = ujj + Auj is handled in the same way. This way, the independent frequency components of the signal 
are set in phase at the focus step by step, i.e. the phase mask e ,A< ^ > ^- ) converges towards the conjugate of e l ^ r ' s ^\ 
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3. Nonlinear detector 

For the first phase, to achieve spatial focusing, the detection of the linear intensity would be sufficient. However, 
for an experimentally convenient approach, we use the nonlinear detector for both spatial and temporal focusing. We 
consider a system with a second-order nonlinearity, which could be realized, e.g., by a nonlinear scatterer, such as 
microbubbles in acoustics or nonlinear nanocrystals [It} in optics, combined with a time-integrated detection of 
the scattered second-order response. We calculate the detector response by 

J SH = aj dt\F rno (t) 2 \ 2 , (26) 

where a is an arbitrary prefactor and F mo (t) the amplitude at the focus point in the time domain. For the narrowband 
pulses emitted in the first phase, we can assume that the linear intensity can be readily obtained by I = V J SH for 
the calculation of the transmission coefficients (Eq. [T3")) . 

III. EXPERIMENT 
A. Transfer matrix measurement 

As a first step of the experiment, we recorded the transfer matrices of a multiple scattering medium over a broad 
range of frequencies. The experimental setup consists of two identical transducers, one used as emitter and the other 
as receiver respectively, placed in water on opposite sides of the medium (Fig. []}. The latter is an arrangement 
(thickness L = 80 mm) of parallel steel rods with density 18 cm -2 and a rod diameter 0.8 mm. Both the emitting and 
the receiving transducer can be translated on a line parallel to the medium, on which we defined each 128 emitter 
points and 128 receiver points spaced by 0.4 mm. From each emitter point, we sent an ultrasonic pulse at a central 
frequency of 3.5 MHz with a relative bandwidth of 31%. In this frequency range the mean free path of the medium is 
I w 4mm as determined in [l8| . The transmitted signals were recorded on all receiver points with a sampling frequency 
of 10 MHz. From the Fourier transform of the temporal signal, we obtain the transfer matrices in the considered 
frequency range (Eq. QJ . We determined the correlation of the transmission coefficients both in the spatial and in the 
frequency domain by 

G,, (Aw) = (T nm (w)T* m (w + Aw)) n)Tn , (27) 

C r (Ar) = (T nm (u>)T: (m+Ar) (u>)) , (28) 

where the average is performed over the denoted matrix coefficients in Eq. [27] and over the first matrix entry and 
frequencies in Eq.[55J From the full width at half maximum of the correlation functions, 5u) = FWHM(|C W |) = 5.14kHz 
and FWHM(|C r |) = 1.04mm, we can determine that we have 1/5. 14kHz independent frequencies per frequency unit 
and N = (128 ■ 0.40mm) /1.04mm = 49.2 independent emitter points available [l||. 

With the transfer matrices at hand, we can calculate the linear response of the system for arbitrary signals emitted 
either from the 'emitter' or 'transducer' side. For all following experiments, we emitted pulses with a Gaussian spectral 
function of 10% bandwidth around the central frequency from all transducers. A typical signal in the receiver plane 
is plotted in Fig. [U(a) and Fig.^a). The energy is spread widely both temporally and spatially. The ballistic part of 
the wave has disappeared, which confirms that the medium is strongly scattering. 

B. Spatial focusing 

Starting from this emission we performed a frequency- resolved wavefront shaping optimization as described in section 
pi B 1[) . The bandwidth of the narrowband pulses is chosen to Aw = 1.4kHz to be smaller than the correlation <5w. 
In two experiments we use both emission with optimal phase and amplitude (Eq. 1161) and the emission with optimal 
phase only (Eq. |2"2")) . The deposed energy in the receiver plane for both experiments is plotted in Fig. @] It is 
calculated by integrating over time the square of the wave amplitude at the receiver points (see Fig. [2jb) for the signal 
obtained in optimizing both amplitude and phase). The optimization leads to a spatial focusing in both cases. For 
phase and amplitude control, the increase in energy at the focal spot is r\ = 49.2, which corresponds perfectly to the 
number of estimated independent emitters N = 49.2 (Eq. I2ip . The enhancement by the time-reversal focus is slightly 
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Figure 2: Received amplitude for nouoptimized emission (a), for spatial wave front shaping (b) and subsequent 
temporal focusing (c), and for comparison, for the time- reversal focus (d). The areas in the dashed rectangles are 

plotted magnified in the inlays. 



higher, indicating remaining deviations between the wave front shaping emission and the time-reversal emission. In 
further calculations, we observed that a further reduction of Aw eliminates these deviations, which shows that the 
single- frequency TR emission and the WFS emission are identical (see Eq. [TH] and Eq. [S]) . 

For phase-only optimization, we observe an enhancement lowered by a factor 0.71 (rj po = 35.1). The factor 0.8 
(Eq. [23)1 is not fully reached, which we attribute to the fact, that due to the geometry of the setup not all emitters 
contribute equally at the focal point, effectively reducing the number of emitter points [HI- The resulting time- 
resolved signals for the first case are shown in Fig. HJb) and Fig. [2Kb). The energy is still spread temporally, since 
temporal focusing is impeded by the remaining random phase relation between independent frequency components 
(Eq. [T4|) . In the next step, we determine and correct for this phase factor. 



C. Spatiotemporal focusing 

Starting with the emission for spatial optimization, we perform frequency phase filtering as described in section 
pi B 2[) . After two cycles over all frequencies, the received signals are not only spatially, but also temporally focused 
(Fig.[2Kc) and Fig.[3lc)). The resulting peak intensity can be put in context with the degrees of freedom available for the 
optimization. We have N spa u a i = 49.2 independent emitters at hand. In the frequency domain, Nf req = 96.3 degrees 
of freedom are available calculated from the frequency correlation, Sui = 5.14kHz and a bandwidth of Ac^f, = 495kHz. 
Hence, the total number of degrees of freedom is rjtheory — N spat i a i • Nf req = 4.7 ■ 10 3 . We estimate the intensity 
enhancement to 7] exp = 6.2 • 10 3 , comparing the peak intensity to the intensity before optimization (obtained from the 
transmission averaged over the transducers around the focus around the peak of the diffuse transmission, smoothed by 
a 60/iS window to remove temporal speckle). This number corresponds well to the order of magnitude of the number 
of degrees of freedom. 

As a reference, we performed an equivalent time-reversal experiment, which intrinsically reaches spatiotemporal 
focusing (Fig.[5](d) and Fig.[3Jd)). The quality of the focusing from our adaptive wave front shaping and time- reversal 
are nearly identical (Fig. [5]). The small deviations are a result of small remaining phase differences which can be 
minimized by a reduction of the frequency steps during the spatial optimization (see above) and further iterations 
of the adaptive phase filtering process. As an alternative to the step-by-step algorithm we used here, so-called 
stimulated annealing algorithms [l5[ or genetic algorithms [2(| should equivalently find the optimal phase relation. 
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Figure 3: Received amplitude at the focal point for nonoptimized emission (a), for spatial wave front shaping (b) 
and subsequent temporal focusing (c), and for comparison, for time-reversal focusing (d). 
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Figure 4: Received energy before (dotted line), after wave front shaping with phase-only control (dash-dotted line) 
and combined control of amplitude and phase (solid line). The focal spot obtained with time-reversal is plotted for 

comparison (dashed line). 



These algorithms are known to be robust in the presence of noise, but have the disadvantage to turn inefficient for an 
increasing number of degrees of freedom. 



IV. CONCLUSIONS 



In conclusion, we have presented a new approach for spatiotemporal focusing through complex scattering media 
by feedback-based wavefront shaping. In contrast to previous works, our approach is capable of reaching the limit of 
a spatiotemporal matched filter; the maximum possible amount of the emitted energy is deposed at the target. We 
showed that to achieve this, phase and amplitude of the emission need to be controlled at each frequency. In contrast 
to the classical time reversal experiment that achieves the ideal focus, here the direct access to the field amplitude 
at the target is not required as only the intensity associated to a nonlinear response needs to be detected. This 
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Figure 5: Amplitude after wave front shaping (solid line) and time reversal (dashed line) on the focal spot. 



point of is particular interest for many applications in wave physics where the accessible information is restricted to 
field intensity. For example, in optics, our method can be realized by fluorescent dye molecules as they are used for 
two-photon microscopy. In ultrasound therapy, it could optimize both spatial and temporal focusing of therapeutic 
beam through the skull bone thanks to MR radiation force imaging. Being generally applicable to all types of waves, 
we believe that our method is promising for a wide range of applications in imaging and sensing, and for the control 
of wave propagation in combination with complex media such as new metamaterials. 
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of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is supported financially by Nederlandse 
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